### Country-wise drop 
rm(list = ls())


data <- read_excel("Data/2_data_final_std.xlsx")


data_europe <- data %>% 
  subset(region=="Europe" &  (polity>=6 | is.na(polity))) #Ukraine drops out

### Country-wise drop 

library(lme4)

data1 <- data_europe%>% 
  subset(!is.na(restrict_freedom_index)&!is.na(DQ_new) & country_DB!="Kosovo")

data2 <- data_europe%>% 
  subset(!is.na(restrict_freedom_index)&!is.na(INDLIB_new) & country_DB!="Kosovo")

data3 <- data_europe%>% 
  subset(!is.na(restrict_freedom_index)&!is.na(RULEOFLAW_new) & country_DB!="Kosovo")

data4 <- data_europe%>% 
  subset(!is.na(restrict_freedom_index)&!is.na(MUTUCONS_new) & country_DB!="Kosovo")

data5 <- data_europe%>% 
  subset(!is.na(power_concentration_index)&!is.na(DQ_new) & country_DB!="Kosovo")

data6 <- data_europe%>% 
  subset(!is.na(power_concentration_index)&!is.na(INDLIB_new) & country_DB!="Kosovo")

data7 <- data_europe%>% 
  subset(!is.na(power_concentration_index)&!is.na(RULEOFLAW_new) & country_DB!="Kosovo")

data8 <- data_europe%>% 
  subset(!is.na(power_concentration_index)&!is.na(MUTUCONS_new) & country_DB!="Kosovo")

countries1 <- unique(data1$country_DB)
countries2 <- unique(data2$country_DB)
countries3 <- unique(data3$country_DB)
countries4 <- unique(data4$country_DB)

countries5 <- unique(data5$country_DB)
countries6 <- unique(data6$country_DB)
countries7 <- unique(data7$country_DB)
countries8 <- unique(data8$country_DB)


coef.catcher1 <- matrix(NA,length(countries1),1)
var.catcher1 <- matrix(NA,length(countries1),1)

coef.catcher2 <- matrix(NA,length(countries2),1)
var.catcher2 <- matrix(NA,length(countries2),1)

coef.catcher3 <- matrix(NA,length(countries3),1)
var.catcher3 <- matrix(NA,length(countries3),1)

coef.catcher4 <- matrix(NA,length(countries4),1)
var.catcher4 <- matrix(NA,length(countries4),1)

coef.catcher5 <- matrix(NA,length(countries5),1)
var.catcher5 <- matrix(NA,length(countries5),1)

coef.catcher6 <- matrix(NA,length(countries6),1)
var.catcher6 <- matrix(NA,length(countries6),1)

coef.catcher7 <- matrix(NA,length(countries7),1)
var.catcher7 <- matrix(NA,length(countries7),1)

coef.catcher8 <- matrix(NA,length(countries8),1)
var.catcher8 <- matrix(NA,length(countries8),1)


for (i in countries1){
  ###  Models
  
  mod1 <- lm(restrict_freedom_index ~ DQ_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher1[which(i==countries1),1] <-  coef(mod1)[2]
  var.catcher1[which(i==countries1),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries2){
  ###  Models
  
  mod1 <- lm(restrict_freedom_index ~ INDLIB_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher2[which(i==countries2),1] <-  coef(mod1)[2]
  var.catcher2[which(i==countries2),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries3){
  ###  Models
  
  mod1 <- lm(restrict_freedom_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher3[which(i==countries3),1] <-  coef(mod1)[2]
  var.catcher3[which(i==countries3),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries4){
  ###  Models
  
  mod1 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher4[which(i==countries4),1] <-  coef(mod1)[2]
  var.catcher4[which(i==countries4),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries5){
  ###  Models
  
  mod1 <- lm(power_concentration_index ~ DQ_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher5[which(i==countries5),1] <-  coef(mod1)[2]
  var.catcher5[which(i==countries5),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries6){
  ###  Models
  
  mod1 <- lm(power_concentration_index ~ INDLIB_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher6[which(i==countries6),1] <-  coef(mod1)[2]
  var.catcher6[which(i==countries6),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries7){
  ###  Models
  
  mod1 <- lm(power_concentration_index ~ RULEOFLAW_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher7[which(i==countries7),1] <-  coef(mod1)[2]
  var.catcher7[which(i==countries7),1] <-  vcov(mod1)[2,2]
  
}

for (i in countries8){
  ###  Models
  
  mod1 <- lm(restrict_freedom_index ~ MUTUCONS_new + hospital_beds + gdp_pc, 
             data =  data_europe[data_europe$country_DB!=i,])
  
  coef.catcher8[which(i==countries8),1] <-  coef(mod1)[2]
  var.catcher8[which(i==countries8),1] <-  vcov(mod1)[2,2]
  
}

z.catcher1 <- coef.catcher1/sqrt(var.catcher1)
z.catcher2 <- coef.catcher2/sqrt(var.catcher2)
z.catcher3 <- coef.catcher3/sqrt(var.catcher3)
z.catcher4 <- coef.catcher4/sqrt(var.catcher4)
z.catcher5 <- coef.catcher5/sqrt(var.catcher5)
z.catcher6 <- coef.catcher6/sqrt(var.catcher6)
z.catcher7 <- coef.catcher7/sqrt(var.catcher7)
z.catcher8 <- coef.catcher8/sqrt(var.catcher8)



for (k in 1:3){
  countries <- eval(parse(text=paste0("countries",k)))
  coef.catcher <- eval(parse(text=paste0("coef.catcher",k)))
  var.catcher <- eval(parse(text=paste0("var.catcher",k)))
  mp <- coef.catcher[,1]-1.96*sqrt(var.catcher)[,1]
  pp <- coef.catcher[,1]+1.96*sqrt(var.catcher)[,1]
  max.pp <- max(pp)
  min.pp <- min(mp)
  span <- max.pp-min.pp
  if (k==1)   pdf("Out/FigureA1_app.pdf",width = 12,height = 6) 
  if (k==2)   pdf("Out/FigureA2_app.pdf", width = 12,height = 6) 
  if (k==3)   pdf("Out/FigureA3_app.pdf",width = 12,height = 6) 
  #pdf(paste("Out/rplot",k,".pdf",sep = "")) 
  par(oma=c(6,1,0,0))
  plot(0, xaxt = 'n',  bty = 'n', pch = '', ylab = '', xlab = '', 
       xlim=c(0,length(countries)), ylim=c(min.pp-0.1*span,max.pp+0.1*span))
  for (i in 1:length(countries)){
    points(i,coef.catcher[i,1], pch=19)
    segments(i,coef.catcher[i,1]-1.96*sqrt(var.catcher)[i,1],i,coef.catcher[i,1]+1.96*sqrt(var.catcher)[i,1], lwd=0.5)
    segments(i,coef.catcher[i,1]-1.64*sqrt(var.catcher)[i,1],i,coef.catcher[i,1]+1.64*sqrt(var.catcher)[i,1], lwd=1.5)
  }
  abline(h=0, lwd=2)
  axis(side =1, at = c(1:length(countries)),las = 3,labels = countries)
  dev.off()
}
